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Abstract We compute the energy dependence of the Pj- 
integrated cross section of directly produced quarkonia in 
pp collisions at next-to-leading order (NLO), namely up 
to a- 3 , within nonrelativistic QCD (NRQCD). Our analy¬ 
sis is based on the idea that the Pj-integrated and the Pj- 
dijferential cross sections can be treated as two different ob¬ 
servables. The colour-octet NRQCD parameters needed to 
predict the Pj -integrated yield can thus be extracted from 
the fits of the Cy -differc ntial cross sections at mid and large 
Pj. For the first time, the total cross section is evaluated 
in NRQCD at full NLO accuracy using the recent NLO fits 
of the Fy -differential yields at RHIC, the Tevatron and the 
LHC. Both the normalisation and the energy dependence of 
the J/ip, t// and T(l.S'), we obtained, are in disagreement 
with the data irrespective of the fit method. The same is true 
if one uses CEM-like colour-octet NRQCD parameters. If, 
on the contrary, one disregards the colour-octet contribution, 
the existing data in the TeV range are well described by the 
a 3 contribution in the colour-singlet model -which, at a 4 s , 
however shows an unphysical energy dependence. A similar 
observation is made for This calls for a full NNLO or 
for a resummation of the initial-state radiation in this chan¬ 
nel. In any case, past claims that colour-octet transitions are 
dominantly responsible for low -Pj quarkonium production 
are not supported by our results. This may impact the inter¬ 
pretation of quarkonium suppression in high-energy proton- 
nucleus and nucleus-nucleus collisions. 


1 Introduction 

Understanding the production mechanism of low-By 
quarkonia in nucleon-nucleon collisions is of fundamental 
importance to properly use them as probes of deconfinement 
or collectivity in heavy ion collisions. Indeed, most of the 
analysis of quarkonium production in nucleus-nucleus colli¬ 


sions are carried out on the bulk of the cross section, namely 
at low Pj. 

Recently, comparisons between ALICE data [1] without 
Pj cut and CMS data [2] with Pj cut in PbPb collisions 
at y,v iV :V = 2.76 TeV showed an unexpected suppression 
pattern of the charmonia, at variance with the simple pic¬ 
ture of quarkonium melting in deconfined quark matter [3], 
However, to properly interpret this observation, it is essen¬ 
tial to rule out the possibility that a part of the effect ob¬ 
served could be due to a difference in the production mech¬ 
anism in individual nucleon-nucleon collisions at low and at 
larger Pj. The propagation of a colour-octet pair in a decon¬ 
fined medium certainly differs from that of a colour-singlet 
pair; this can result into a different nuclear suppression (see 
e.g. [6]). On the contrary, as regards the bottomonia, the ob¬ 
servation of the expected sequential-suppression pattern has 
been claimed by CMS [4, 5]. 

Further, the effect of normal nuclear matter may also sig¬ 
nificantly depend on how the pair is produced: the recently 
revived fractional energy loss [7, 8] would for instance act 
on long-lived colour-octet states and probably differently if 
the heavy quark state is already produced colourless at short 
distance, as postulated in the CSM [9]. Saturation effects in 
pA collisions also do depend on the colour state of the per- 
turbatively produced heavy-quark pair [10-12] 

Despite the possibility that NRQCD factorisation would 
not hold at low Pj , several NRQCD analyses have thus 
been carried earlier to evaluate the impact of the colour- 
octet channels to the Pj -integrated J /(// yields [13-15]. A 
first study of the impact of initial state radiations (ISR) on 
the very low Pj J/i/ds and T’s was recently carried out suc¬ 
cessfully in NRQCD [16] - yet at the cost of introducing 
additional non-perturbative parameters. 

Whereas, based on an analysis of the sole early RHIC 
data. Cooper et al. argued [14] that the universality of 
NRQCD was safe and that colour-singlet contributions to 
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the Pj -integrated J/ip yields were negligible, the global 
analysis of Maltoni et al. at NLO showed [15] that the 
colour-octet Long-Distance Matrix Elements (LDMEs) re¬ 
quired to describe the total prompt 7/t/r yield from fixed- 
target energies to RHIC were one tenth of that expected from 
the - leading-order- fit of the R/ -differential cross sections 
at Tevatron energies. 

Such fits of the ^/ -differential J/i/z cross sections have 
recently been extended to NLO —i.e. one-loop- accuracy on 
the prompt 7/t // yields -some of them focusing on the larger 
Pj data and explicitly including the feed-down contribu¬ 
tions [17, 18], some enlarging the analysis beyond hadro- 
production and including rather low-TV data [19]- and on 
the T(nS) yields [20, 21], Thanks to these studies, we can 
significantly extend the existing NRQCD studies of the Re¬ 
integrated cross section by combining in a coherent man¬ 
ner, the hard parts - or Wilson coefficients- up to a $, first 
derived by [22] and which we have systematically checked 
with FDC [23], with the NRQCD matrix elements fitted at 
NLO on the Pj dependence of the yields. One can indeed 
consider the P j--integrated and the Pj-differential cross sec¬ 
tions as two different observables - their Born contributions 
involve different diagrams - and such a procedure is not a 
all trivial physics-wise. 

As we detail later, our results show that the data do 
not allow for a global description of both the A/--integrated 
and R/ -differential quarkonium yields. As a point of com¬ 
parison, we also had a look at Colour-Evaporation-Model- 
like (CEM) predictions derived from NRQCD following the 
work of [26] and we found out that it cannot reproduce Re¬ 
integrated yields using the LDMEs obtained following the 
relations of [26] after identifying the minimal singlet transi¬ 
tion to that of the CSM. A contrario , results obtained from 
the traditional CEM implementation at one loop do not show 
a similar issue. 

The inability of colour-octet dominance within NRQCD 
to provide a global description of both low and large Re data 
is in line with the recent findings [28-31] that the sole LO 
colour-singlet contributions are sufficient to account for the 
magnitude of the total cross section and its dependence in 
rapidity, d<r/dy, from RHIC, Tevatron all the way to LHC 
energies. Any additional contribution in this energy range 
creates a surplus 1 as compared to data. 

However, as we also study in a dedicated section, the 
total NLO CSM cross section shows a weird energy depen¬ 
dence at LHC energies. The problem is striking for the J/i/z, 
less for the T. In any case, one should be very careful in in¬ 
terpreting these results. In particular, such NLO results can¬ 
not be considered as a improvement of the LO ones. We also 
observe the same issue for r/ c and pb production for which 
there is no final-state-gluon radiation at Born order. We are 

*We however note that the CO contributions by themselves can even 
also overshoot the data. 


therefore tempted to attribute this behaviour to large loop 
contributions which become negative at low Pj, rather than 
to specific effects related to the \S \ production per se. A 
quick inspection of the rather concise one-loop results [32] 
for Tj c and ///, production in the TMD factorisation formalism 
unfortunately does not reveal any obvious negative contri¬ 
butions and does not help in the understanding of this rather 
general issue of quarkonium production in collinear factori¬ 
sation. 

The same problems appear with some CO channels as 
well and may therefore cast doubts on the reliability of our 
results in the ffs region where some contributions shows 
a strange behaviour -in particular at LHC energies. At this 
stage, we are not able to conclude from our observations 
whether these problems are indicative of the break down of 
NRQCD factorisation at low Pj and low x or at low x only. 
However, for sure, none of the above observations can rea¬ 
sonably support the idea that CO transitions are dominant at 
low Pj. Such a conclusion would be for the least premature. 

The structure of the paper is as follows. In section 2, we 
detail the procedure to evaluate the Pj -integrated yield at 
one-loop accuracy in NRQCD and we explain the idea un¬ 
derlying this first complete one-loop analysis. In section 3, 
we explain our selection of LDMEs determined at NLO. In 
section 4, we briefly comment on the existing world data sets 
for J/t/z, iff2S) and T(IS) 2 . We also explain how we esti¬ 
mate the direct yields. In section 5, we show and discuss our 
results for the first full one-loop NRQCD analysis of quark¬ 
onium hadroproduction. To go further in the interpretation 
of some of our results, we discuss in section 6 the prediction 
of NRQCD using CEM-like LDMEs. This is also compared 
with the conventional approach based on quark-hadron dual¬ 
ity. Section 7 focuses on CSM results both for the \S \ states 
considered here and for the 1 S o states for which analytical 
results exist. Our conclusions are presented in section 8. 

2 A full one-loop cross-section computation 

2.1 Generalities 

Following the NRQCD factorisation, the cross-section for 
quarkonium hadroproduction can be expressed from the par- 
ton densities in the colliding hadrons, f(x ), a hard-part - 
the partonic cross section- for the production of a heavy- 
quark pair with zero relative velocity, v, in a definite angular- 
momentum, spin and colour state, and a LDME connected 

^Whereas the Y analysis of Gong et al. [21] treats the T(IS), Y(2S) 
and T(3S), the lack of knowledge on the^ , /,(2R) andxbffP) yields and 
their corresponding feed down to Y(2S) and T(3S) makes the analy¬ 
sis of their direct yield delicate; we have thus decided not to consider 
these in the present study. Our choice has been confirmed by the recent 
LHCb result [63] that a large fraction of the Y(2S) and T(3S) yield ac¬ 
tually comes from / y/,(2R) and XbOP) decays - up to 40% in the T(35) 
case. 
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Fig. 1: Representative diagrams contributing (a-b) at Born order to i + j —> <3, (c-e) both at Born order to i + j —> <3+jet and at 
one loop to i + j —* <3, (f) at one loop to i + j —> Q, (g-k) at one loop to i + j —> <3+jet. 


to the hadronisation probabilities of the intermediate state 
into the quarkonium. Namely, one has for the production of 
a quarkonium <3 along with some unidentified set of parti¬ 
cles X, 

cr =Yj f dx{dx 2 fifjd-[i + j (QQ)„ + X]{0 Q (n)) (i) 

Uj,n 

where the indices i, j run over all partonic species and n de¬ 
notes the colour, spin and angular momentum states of the 
intermediate QQ pair. 

For the 3 S i quarkonium states, the first CO states which 
appear in the v expansion are the j^, 3 S ^ and 3 P^ states, 
in addition to the leading v contribution 3 .S' ^ from a CS tran¬ 
sition. One however has to note that, for hadroproduction, 
whereas the CO contributions already appear at a ^ (Fig. la 
& lb), the CS one only appear at <x 3 (Fig. lc). These a 3 CO 
graphs nevertheless do not contribute to the production of 
quarkonia with a nonzero Pj since they would be produced 
alone without any other hard particle to recoil on. 

The Born contributions from CS and CO transitions are 
indeed different in nature: the former is the production of 
a quarkonium in association with a recoiling gluon, which 
could form a jet, while the latter is the production of a quark¬ 
onium essentially alone at low Pj. 

Let us now have a look at the « 3 CO contributions 
(Fig. 1 d- If) which are then NLO -or one loop- corrections 
to quarkonium production and which are potentially plagued 
by the typical divergences of radiative corrections. Yet, the 
real-emission a- 3 corrections to CO contributions (Fig. Id & 
le) can also be seen as Born-order contributions to the pro¬ 
duction of a quarkonium + a jet -or, equally speaking, of a 
quarkonium with Pj » Aqcd- As such, they do not show 
any soft divergences for Pj 4 0. These are supposed to be 
the leading contribution to the ^/ -differential cross section 
in most of the data set taken at hadron collider (Tevatron, 
RHIC and LHC). These are now known up to one-loop ac¬ 
curacy, namely up to a ^ (see e.g. [17-21, 34, 35]) (Fig. lg 
& lh). 

It is important to note that one cannot avoid dealing 
with the divergences appearing at o- 3 if one study the Re¬ 
integrated cross section. 


2.2 Different contributions up to rr 3 

At o' 3 , the CO partonic processes are: 

q + q~>QQ?Sf] 

g+g^QQ[ l sf?pf = 0i2 ] 

where q denotes u, d, s. 

At o' 3 , the QCD corrections to the aforementioned chan¬ 
nels include real (Fig. Id & le) and virtual (Fig. If) cor¬ 
rections. One encounters UV, IR and Coulomb singulari¬ 
ties in the calculation of the virtual corrections. The UV- 
divergences from the self-energy and triangle diagrams are 
removed by the renormalisation procedure. Since we follow 
the same lines as [18, 35] where all the procedure is de¬ 
scribed, we do not repeat its description. As regards the real- 
emission corrections, they arise from 3 kinds of processes 
(not all drawn): 

£ + £ -» ‘S O 81 J 81 > 3/5 /J 0 , 2 ] + 8 > 

g + q(q ) -» ee[ 1 5™, 3 5f ] , 3 ^ 8 i a2 ]+^@, (3) 

q + q^QQl^^Mf^+g. 

As usual, the phase-space integrations generate IR singu¬ 
larities, which are either soft or collinear and can be con¬ 
veniently isolated by slicing the phase space into different 
regions. Here we adopt the two-cutoff phase space slicing 
method to deal with the problem [36]. 

As we previously alluded to, the a 3 CS contribution is 
particular since, in the limit v = 0, it would be strictly speak¬ 
ing Born order for both the production of a quarkonium and 
of a quarkonium + a jet. It arises from the well-known pro¬ 
cess: 

g + g^QQ[ 3 S [ "]+g (Fig. lc) (4) 

Our calculations is equivalent to a previous work by 
Maltoni et al. [15, 22] and we have checked that we repro¬ 
duce their results for all the relevant channels. As announced 
in the introduction, one of the novelty in our study resides in 
the use of the LDMEs fitted at the same order, Le. one loop, 
to the Ry -differcntial cross sections. As such, this is the first 
global NLO analysis of hadroproduction. 


(Fig. la) 
(Fig. lb) 
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Since we also look at data at rather low energies, we also 
included a CS channel via y* exchange. Indeed, as noted in 
a different context in [37], the QED CS contributions via 
y* are naturally as large as the CO ] S ^ transition via g* 
- the a em suppression being compensated by the small rel¬ 
ative size of the 3 S ® CO LDME (<9(10 -3 )) as compared to 
the 3 5 ^ CS LDME (0(1 )). The real-emission contributions 
arise from 

q + q->QQ?S { "]+g, 

g + q(q)^QQ[ 3 S\ l] ] + q(q), 

whereas the loop contributions are only from 

q + q-^ QQ[ 3 S [ "l ( 6 ) 

Fig. la (Fig. If) with the s-channel gluon replaced by a y* 
would depict the Born (a one-loop) contribution. We have 
however found that they do not matter in the regions which 
we considered. 

3 Constraints on the LDMEs from the Pj -differential 
cross section 

The CS LDMEs can either be extracted from the leptonic 
decay width at NLO or can be estimated by using a poten¬ 
tial model result, which gives for the Buchmuller-Tye po¬ 
tential [38] \R JH ,m 2 = 0.81 GeV 3 , |^ (2S) (0)| 2 = 0.53 GeV 3 
and |flrus)(0)| 2 = 6.5 GeV 3 . 

As regards the CO LDMEs, they can only be extracted 
from data. As we discussed above, our aim is to analyse the 
Pj-integrated yield using the constraints from the Pj de¬ 
pendence of the yields. 

3.1 J/tp 

In the J / if/ case, we will use the results of five fits of this de¬ 
pendence [17, 18, 39-41]. The first two were limited to pp 
data but explicitly took into account the effect of the feed- 
down 1 The latter fit was based on a wider set of data includ¬ 
ing ep and yy systems but the feed-down effects were only 
implicitly included through constant fractions for these sys¬ 
tems. The fourth one includes the recent r/ c measurement at 
Pt > 6 GeV by LHCb [42] by relying on the heavy-quark 
spin symmetry of NRQCD which relates colour-octet ma¬ 
trix elements of spin-singlet and triplet quarkonia with the 
same principal quantum number n. The fifth one incorpo¬ 
rates the leading-power fragmentation corrections together 
with the usual NLO corrections, which results in a different 
short-distance coefficient and allows for different LDMEs. 

3 In [17], Ma et al. used both the prompt J/xj/ yield and polarisation 
data from CDF(run II). In [18], Gong et al. chose to fit the CDF and 
LHCb experimental data for the yield only (no polarisation data). 


Table 1: Values of (Oj^( 3 P [3] )), and 

(0.,/^S ' X ')) from 5 NLO ( i.e. at one loop) fits of the 
Pt dependence of the yields, which we will use to compute 
[the energy dependence of] the Pt -integrated yields. 


Ref. 

(O m ( 3 pf)) 

(in GeV 5 ) 

(O/aA'sI, 81 )) 

(in GeV 3 ) 

(Oj M ( 3 s[ 8] )) 

(in GeV 3 ) 

[17] 

2.1 x 10 -2 

3.5 x 10 -2 

5.8 x 10 -3 

[40] 

3.8 x 10 -2 

0.7 x 10“ 2 

1.0 x 10 -2 


3.4 xlO' 2 

0.9 x 10 -2 

1.6 x 10 -2 


4.3 xlO- 2 

0 

1.1 xl0‘ 2 


4.5 x 10 -2 

1.6 x 10 -2 

1.2 xlO' 2 


5.4 x 10 -2 

0 

1.4 xlO' 2 


2.3 x 10 -2 

1.6 x 10 -2 

0.6 xlO' 2 


3.2 x 10 -2 

0 

0.8 xlO' 2 

[18] 

-2.2x 10 -2 

9.7 x 10“ 2 

-4.6 x 10 -3 

[39] 

-9.1 x 10 -2 

3.0 x 10 -2 

1.7 x 10 -3 

[41] 

1.1 xlO' 2 

9.9 x 10 -2 

1.1 xlO' 2 


Another recent fit [43] took the t] c measurement into ac¬ 
count. The LDME values which they found fall into range 
considered for [40], therefore we do not use it separately. 

In Ref. [17], Ma et al. have have based their analyses on 
the fit of two linear combinations 4 of LDMEs: 

+ ^{Oj^pf)), 

nic (7) 

M'W - (O m ( 3 S [ , 8] )> + -^{Oj^Pf)). 

m c 

They proceeded to two fits with different Pt cuts. We use 
that for Pt >7 GeV and limit ourselves to the central val¬ 
ues they obtained: = 7.4 x 10 -2 GeV 3 and = 

J o ,r 0 l,n 

0.05 x 10 -2 GeV 3 , since a single set of values of and 

translates anyhow into a wide range of values of the 
LDMEs. Indeed, limiting ourselves to positive values of 
(Oj^Sf)) and (Oji^S j 8 ^)), one can solve Eq. 7 and get 
the loose constraint: (Oj^( 3 P [3] )) e [-0.2,4.1] x 10“ 2 GeV 5 . 
As a central value, we choose the middle of the allowed 
interval. The same group has however recently improved 
their analysis by taking into account the feed-down [44]. As 
aforementioned, they in turn performed a new fit [40] in¬ 
cluding ri c data. The 6 sets of LDMEs to be used to probe 
the allowable parameter space of the fit are given in Table 1. 

As mentioned above, in [39], Butenschoen et al. pro¬ 
ceeded to a global fit of prompt J /i/y data from pp, yy, yp 
systems 5 . Since yy, yp mostly lies at low Pt, they also 
considered data at rather low Pt from RHIC. They did 
not included NLO predictions for Xc in the fit. Rather they 
assumed a constant direct fraction, for instance 36 % for 
hadroproduction. 

4 r 0 =4.1 and n =-0.56 at 7 TeV 

5 and one point from e + e~ at KEKB. 
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3.2 MIS) 

Buttenschoen et al. did not provide a fit of 1^(25) in [39] 
due to the lack of data besides those from pp collisions. 
The LDMEs which we consider for i//(2S ) are therefore 
only from [17] and [18]. For the former fit, the values are 
obtained in the same way as for the J/tp, where = 

2.0 x 10“ 2 GeV 3 and = 0.12x 10“ 2 GeV 3 . The result¬ 

ing values as well as those from [18] are gathered in Table 2 


Table 2: Same as Table 1 for t/425). 


Ref. 

(^( 2 S)( 3 T™)> 
(in GeV 5 ) 

(cw 1 ^ 81 )) 

(in GeV 3 ) 

<a«2.S)(M 8] )> 

(in GeV 3 ) 

[18] 

9.5 x 10 -3 

— 1 . 2 x 10~ 4 

3.4 x 10 “ 3 

[17] 

-4.8 xlO ^ 3 

2.9 x 10 -2 

0 


7.9 x 10 -3 

5.6 xlO ' 3 

3.2 x 10 “ 3 


1.1 xl 0~ 2 

0 

3.9 x 10 “ 3 


In [44], the authors of [17] tried to refit the existing data 
with a larger Pj cut-ofF. Such a fit already badly overshoots 
mid-/ 7 / data. We therefore do not consider it in this work. 
For the same reason, we have not considered the fit of [45] 
since it only reproduces the 1/425) data in an admittedly nar¬ 
row -high Pj- range. 


3.3 T(15) 

As regards the T(15), there are two NFO analyses from [20] 
and [21], However, Wang et al. used in [20] a different 
value of the NRQCD factorisation scale p.\ which we use 
in the present evaluation, that is p A — mi,. To perform a cor¬ 
rect comparison would have required a new evaluation of the 
hard coefficients with their choice of p A to use their FDME 
values. In addition, although they did consider the effects 
of excited feed-down, they have not disentangled the direct 
contribution to that of the feed-down in their FDME extrac¬ 
tion. The central values of [21] are gathered in Table 3. 


Table 3: Same as Table 1 for T(15). 



Ref. 

<0r(is)( 3 T[, 8] )> 

(Oras/Sl, 81 )) 

(Dr<ts)( 3 Sj 8] )> 


(in GeV 5 ) 

(in GeV 3 ) 

(in GeV 3 ) 

[ 21 ] 

-13.6x 10 -2 

11 . 2 x 10 -2 

-4.1 x 10 -3 



4 World data and feed-down effects 

As regards the data for J /<// and 14(25), we drew on the ex¬ 
tensive set used in [15] with the exception that we only kept 
data: 

- derived from more than 100 events at a given Vv 

- from pp or pp collisions only in order to avoid dealing 
with nuclear effects; 

- where dcr/dy was derived at y — 0. 

Table 4: 7/t4 data set used in our data-theory comparison. 
The experimental values quoted in the experimental papers 
may have been multiplied by feed-down factors (see text). 


Experiment/Collaboration 

M (GeV) 

^^-extr.direct 

dy 

n (nb) 

y=0 

UA 6 [46] 

24.3 

3.7: 

+1 

ISR-Clark et al. [47] 

52.4 

6.7 ± 

1.2 

ISR-R806 [48] 

53 

8.2 ± 

2.3 

ISR-Clark et al. [47] 

62.7 

6.2 ± 

1.2 

ISR-R806 [48] 

63 

9.0 ± 

2.5 

PHENIX [49] 

200 

27 + 6 

CDF [52] 

1960 

27 + 

16 

20 

ALICE [50] 

2760 

138±46 

ALICE [51] 

7000 

220 ±53 


To this set, we have added data published later than 2006 
which includes data from the FHC. We have also added one 
point from the CDF collaboration at the Tevatron' 1 . All the 
quoted uncertainties are combined in quadrature together 
with that of the direct fraction 7 which we assumed to be 
energy independent and FjJ^ ct = 60+ 10% [28]. 

As regards the i4(25), the data sets are very scarce es¬ 
pecially if one focuses on Pj -integrated yields at y = 0. In 
fact, there is only data from ISR-Clark et al. [47] averaged 
over Ms = 52.4 and 62.7 GeV and from PHENIX at \fs - 
200 GeV. CDF measured the cross section at -\As = 1 .96 TeV 
for [y[ < 0.6 but only for Pj > 2 GeV [53]. In order to use 
this precise measurement, we have extrapolated it by assum¬ 
ing the same ratio {y=o/ *^_> 2 GeV) ^ = 0 82 as 

for the JI ijj [52], As for now, there does not exist measure¬ 
ment at FHC energies in the central rapidities down to small 

6 To be precise, the CDF measurements of prompt 7/ i/j did not extend 
lower than Pj = 1.5 GeV, only the sum of prompt and non-prompt 
7/i/f was measured down to Pj = 0. In order to derived a Pj -integrated 
prompt yield, we have made the reasonable hypothesis that the prompt 
fraction was similar below Pj = 1.5 GeV than just above. This induces 
an uncertainty which is certainly irrelevant for the present comparison. 

7 For the LHC data at low Pj, in particular the ALICE data, 90% of 
the yield is considered to be prompt. For all the other measurements - 
mainly at low energies- which did not separate out the prompt and non¬ 
prompt, we assumed the fraction of non-prompt 7/i/r to be negligible 
given the other uncertainties. 
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Fig. 2: (Colour online) The cross section for direct (a) J/i/z, (b) if/(2S) and (c) T(1.S’) as a function of \/V. The blue dot-dashed 
curve is the central CS curve. Its relative uncertainty is shown in the lower panels; the light green (light blue) band shows the 
scale (mass) uncertainty. The dashed red curve is the total CO contribution from 3 channels: 3 / J J ) X ' (thin dot-dashed orange), 
'.S' (thin dotted magenta) and IS ® (thin dashed green). The total CO uncertainty relative to the CS central curve is shown in 
the lower panels; the light pink (purple) band shows the scale (mass) uncertainty. The black is the total contribution (CS+CO) 
at one loop. These are compared to experimental data (see text) multiplied by a direct fraction factor (when applicable) and 
normalised to the central CS curve in the lower panels. [Negative CO contributions are indicated by arrows]. 
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enough Pj to perform a model-independent enough extrap¬ 
olation 8 . 


Table 5: if/(2S) data set used in our data-theory comparison. 



Experiment/Collaboration 

Vi (GeV) 


ISR-Clark et al. [47] 
PHENIX [49] 

CDF [53] 

52.4-62.7 

200 

1960 

0.2 ±0.07 

0.91+0.23 

4.0 ±0.5 



Table 6: T( 15) data set used in our data-theory comparison. 
[A star indicate that the measurement could not resolve the 
15,25 and 35 states.] 


Experiment/Collaboration yfs (GeV) 


j^extr.direct 

dy y=0 


(pb) 


E866 [54] 

38 

1 . 1 + 0.1 

ISR-R806 [48] 

63 

9.1 ± 3.6 

*STAR [55] 

200 

38.4+ 12.4 

* UA1 [56] 

630 

120 + 36 

CDF [57] 

1800 

380 ±60 

DO [58] 

1960 

410+80 

CMS [2] 

2760 

610+ 170 

ATLAS [59] 

7000 

1180 + 200 

CMS [60] 

7000 

1330 + 230 


As regards the Y(15), the data set is surprisingly wider 
than that of if/(2S) despite a significantly smaller produc¬ 
tion cross section. It is certainly due to the larger energy 
of the decay leptons and to the smaller background. For a 
long time, it was considered that only half of the (low-/ J yj 
T(15) were directly produced =50+ 10%) based 

on an early CDF measurement [61]. Recent LHCb studies 
of Xb production [62, 63] along with T(25,35) cross sec¬ 
tion measurements [59, 60, 64, 65], rather indicate that two 
thirds of the T(15) are directly produced, we will therefore 
opt for = 66 ± 10%. Yet, a number of experiments 

could not resolve the 3 T states. In this case, one should 
apply [28] a slightly smaller direct fraction which we take 
to be Fyvs+2S+3S) =60+ 10%. As we take this fraction to 
be energy independent, we chose a conservative estimate of 
their uncertainty. 

Table 4 shows the J/i/z data set. Table 5 that of i//(2S ) and 
Table 6 that of T. 


8 One could however use the LHCb and ALICE measurements in the 
forward region since the rapidity dependence is certainly is better con¬ 
trol than the Pj dependence from 6 GeV downwards. 


5 Complete NLO results within NRQCD 

In the numerical computation at NLO, the CTEQ6M 
PDF [66] 9 , and the corresponding two-loop QCD coupling 
constant a s are used 10 . The charm quark mass, m c , is set by 
default to 1.5 GeV and the bottom quark one, mb, to 4.5 GeV. 
Our default choices for the renormalisation, factorisation, 
and NRQCD scales are pr- pf = Po with po = 2m q and 
Pa = 111 Q, respectively. When other choices are made, in par¬ 
ticular to estimate the theoretical uncertainty due to the lack 
of knowledge of corrections beyond NLO, they are indicated 
on the corresponding plots. We have taken 5 S = 10 -3 and 
6 C = As /50 for the two phase space cutoffs -the insensitivity 
of the result on the chosen values for these cut-off has been 
checked. Our results for direct J/if/, if/(2S) and Y(15) are 
shown on respectively Fig. 2 (a), (b) and (c). 

We first discuss the comparison between the five fits and 
the J/ip data (Fig. 2 (a)). Without a surprise, our study shows 
that the global fits including rather low Pj I niQ data, that is 
the one of Butenschoen et al. [39] provides the only accept¬ 
able description of the total cross section. We however note 
that the latter fit does not provide a good description of the 
JI if/ polarisation data and, as recently noted [68], it yields 
to negative cross section for Jfifi + y at large Pj. Finally, 
it does not allow [69] to describe the r/ c data. The fits of 
Gong et al. [18], and Ma et al. [17, 40] greatly overshoot the 
data in the energy range between RHIC and the Tevatron, 
whereas these fits a priori provide a good description of the 
/ J / -differcntial cross section at these energies. 

The fit of Bodwin et al. [41] gives the worse account 
of the Pj -integrated Jf ij/ data in the whole energy range. In¬ 
deed, the new ingredient of [41] allows one to describe high- 
P T data with a large 1 5[f 1 LDME (see Table 1) -as for [18] 
but without negative LDMEs for the other octet LDMEs- 
which results in too large a yield at low Pj. 

In addition, we also note the strange energy dependence 
of at least the P -wave octet channel. We postpone its dis¬ 
cussion to section 5.1 where this is analysed in more detail 
for the 1 5 J. 8 ' transition and, in section 7, where we discuss a 
similar observation for the CSM yield at NLO. 

As regards the if/{2S) (Fig. 2 (b)), our NLO NRQCD re¬ 
sults do not reproduce the data at all at RHIC energies and, 
since both fits as dominated by the P -wave octet channel, 
shows a nearly unphysical behaviour at LHC energies. 


9 We have checked by using MSTW [67] that our results do not quali¬ 
tatively change when another PDF set is used. 

10 For the channels which are only considered at tree/Bom level, we 
used the LO PDF set CTEQ6L and the coupling at one loop. Such a 
choice is a matter of convention since no divergence is cancelled be¬ 
tween this contribution and the other contributions. One could have 
chosen a NLO PDF set and a s at two loop. This remark does not con¬ 
cern the real-emission radiative corrections of a given channel which 
are treated as usual, i.e. with NLO PDFs. 



















The comparison for the T(1S) (Fig. 2 (c)), is more en¬ 
couraging. At RHIC energies and below, the agreement is 
even quite good, while at Tevatron and LHC energies, the 
NLO NRQCD curves only overshoot the data by a factor of 

2 . 

We finally note that from RHIC to LHC energies, the 
LO CSM contributions (the blue in all the plots) accounts 
well for the data. The agreement is a bit less good for if/(2S ) 
if we stick only to the default/central value. This is not at 
all a surprise and is in line with the previous conclusions 
made in [28-31], In fact, strangely enough, it seems that 
it is only at low energies (below \[s - 100 GeV) that the 
CO contributions would be needed to describe the data. The 
more recent data from the LHC and the Tevatron tend to 
agree more with the LO CSM. 

Overall, this shows -unless the resummation of ISR 
modifies our predictions by a factor of ten- that it would 
be difficult to achieve a global description of the total and 
TV-differential yield and its polarisation at least for the char- 
monia. 

As we discussed in the introduction, a first resumma¬ 
tion study has recently been performed within NRQCD [16]. 
When combined with the results of [17], this resummation 
yields [16] to a good description of the low-TV data. It 
should however be stressed that this study introduces 3 new 
parameters g 12,3 to parametrise the so-called W NP function 
used the CSS resummation procedure. Moreover, we stress 
that such values of the CO LDMEs would result in a neg¬ 
ative NLO /^ -differential cross section for J/i/z + y at large 
Pj where NRQCD factorisation should normally hold. 

To close the discussion of the theory-data comparison, 
let us note that the 3 S ^ channel alone would provide a de¬ 
cent energy dependence. If we were to refit the low-TV data 
and thus obtain a dominance of the s S j 8 * channel, the yield 
at large Pj would nevertheless dominantly be transversely 
polarised in disagreement with existing data ( e.g. [70]). Yet, 
the better energy dependence of IS' ® at NLO with the re¬ 
spect to the other octet channels, which shows a flat energy 
dependence in the TeV region, is probably a fortunate “acci¬ 
dent”. Indeed, most of the \S ^ yield up aj is in fact not at 
one loop, but for the -suppressed- qq contribution, since the 
gg —> 3 Sj 8 ^ is zero. The curve for 3 S' | 8 ' - as well as 3 S' ^ 1 ^ - 
shown on Fig. 2 (a) is effectively a Born-order one. Clearly, 
these behave better than the channels where the loop contri¬ 
butions are allowed. 


5.1 Behaviour at high \[s and scale dependence of the 
g 8 ^ contribution 

In view of the observations just made in the previous para¬ 
graph, we have found it useful to analyse more carefully the 
behaviour of one specific channel. We have decided to look 


more carefully at the simplest one, that is from the 1 S { 8 ' 
transition in particular at how the scale choices influence the 
behaviour of the yield at large \[s. Analytical results for the 
hard-scattering partonic amplitude squared can be found in 
appendix C of [22]. We have used them in a small numerical 
code to convolve it with PDFs and checked them with the 
results of FDC. The advantage of using FDC is that we can 
easily cut on Pj and y. 



Fig. 3: The cross section for the production of a J/if/ from 
only a colour octet *5^ cc state as a function of the ems- 
energy for various choice of the mass and scales. 

As one could have anticipated from the band of the 
lower panel of Fig. 2 (a), one observes that for a wide range 
of scale choices, the energy dependence remains extremely 
flat. For some of the choice where ftp > /j, R , one even sees 
that the cross section clearly decreases and becomes neg¬ 
ative - when the yield becomes negative the curve stops. 
This is of course not satisfactory. At this stage, we are not 
able to conclude from our observations -normalisation and 
high energy issue- whether these point at the break down of 
NRQCD factorisation, NRQCD universality or should force 
us to continue questioning our understanding of the mid- 
and high-Pj quarkonium production mechanisms. To inves¬ 
tigate this a bit further, we look at the predictions of another 
approach -the colour evaporation model- in the next sec¬ 
tion and later in the colour-singlet model both for 3 .S' 1 and 
: 5o quarkonia. 


6 Colour-Evaporation-Model-like NRQCD evaluation 

To go further in our investigations of QCD one-loop effects 
on the energy dependence of quarkonium production, we 
have found it useful to compare our results with those of 
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the Colour-Evaporation Model (CEM) which directly fol¬ 
lows from the quark-hadron duality [71, 72], The quark- 
onium production cross section is obtained by considering 
the cross section to produce a QQ pair in an invariant mass 
region compatible with its hadronisation into a quarkonium, 
namely between 2m q and the threshold to produce open 
heavy flavour hadrons, 2 mu- To this, one should multiply a 
phenomenological factor accounting for the probability that 
the pair eventually hadronises into a given quarkonium state. 
Overall, one considers 


(N)LO, direct 

(T q 


= V\ 


direct 

a 


„ , (N)LO 

2m « dCr QQ 

dm Q Q 


n2m H 

J2mQ 


dm Q Q 


( 8 ) 


In a sense, the factor 7 > ^ lect , i.e. the probability (or frac¬ 
tion) of QQ pair in the relevant invariant mass region to di¬ 
rectly hadronise into <3, plays a similar role as the LDMEs 
in NRQCD, except that its size can be guessed. Indeed, it 
is expected [73] that one ninth -one colour singlet QQ con¬ 
figuration out of 9 possible- of the open charm cross sec¬ 
tion in this invariant mass region eventually hadronise into 
a “stable” quarkonium. Taking into account this factor 9, in 
the case of J jif/, it was argued [73] that a simple statistical 
counting, which would give: 

^direct = I 2 ~V + 1 = _L f9) 

9 2,-(2/i + l) 45’ 


where the sum over i runs over all the charmonium states 
below the DD threshold, could describe the existing data in 
the late 90’s. The solid turquoise curve -computed at NLO 
as opposed to the analysis of [73]- on Fig. 4 (a) illustrates 
this n . 

Following the fit of Vogt in [75], !P^ ct lies between 
1.5 % and 2.5 %. This indeed remarkably coincides with 
the simple statistical counting 12 . For the T, the correspond¬ 
ing quantity is of similar size, between 2 % and 5 %, al¬ 
though following the state-counting argument, one may ex¬ 
pect a smaller number than for J/if/. Fet us nevertheless 
stress that a violation of Eq. 9 cannot be used to invalidate 
the CEM since this relation completely ignores phase-space 
constraints. What CEM predicts is that !P)!j rec1 is process in¬ 
dependent. 

In [26], Bodwin et al. studied the connexion between 
the CEM and NRQCD. Following [26] up to v 2 corrections, 
only 4 intermediate QQ states contribute to IS’ i quarkonium 
production in a CEM-like implementation of NRQCD. One 


"it has been obtained with MCFM [74] with m c =1.5 GeV, /Jjf = Hf = 
2 m c and mn = mo- 

12 As discussed in [24, 25], further counting rules involving the P- 
waves do not work and illustrate the limit of the model. 


indeed has: 

<a ls .,( 3 5'f n )> =3 x<03 S| ( 1 s! ) 1] )>, 

(0 3Sl ( I 5[ ) 8] )>=^x<03 S| ( I 5t I] )), ( 10 ) 

(Osj/S™)) =4x<C>3 Sl ( 1 s! ) 1] )>. 

All these nonvanishing LDMEs are then fixed if one 
makes the reasonable assumption that (0\ s f ( \S ^ )) is indeed 
the usual CS LDME, i.e. (2 J + 1) |R(0)| 2 . As compared 
to the results presented in the previous section, the only ad¬ 
ditional piece to perform a full one-loop analysis is the hard 
part for '.S’ J, 1 ^ which normally needs not to be considered for 
3 S 1 production at this level of accuracy in v. Computing it 
with FDC [23] does not present any difficulty. 
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Fig. 4: The cross section for direct (a) 7/t/r and (b) T(15) as 
a function of \fs from NLO NRQCD using the CEM-like 
constrained LDMEs assuming a minimal singlet transition. 
It is compared to the existing experimental measurements 
(see text). 
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Fig. 4 (a) & (b) show the resulting cross sections of 7/^ 
and TllS) production for the relevant channels and their 
sum, to be compared to the world data set used in the pre¬ 
vious section. By construction, the 3 S^ curve is the same 
as in the previous section. One notes that, as for the P-wave 
octet, the '.S'J, 1 ' curve strangely flattens out in the 7/t/f case at 
high energies. We will come back to this in the next section. 

The total CEM-like contribution greatly overshoots the 
data, by a factor as large as 100. This was to be expected 
since (i), following Eq. 10, all the LDMEs are roughly of 
the same size, (ii) the 3 S -11 ' is roughly compatible with the 
data and (iii) the hard part for the other transitions appear at 
a j, are not suppressed as Pj —> 0 and are thus expected to 
give a larger contribution than the IS’ ^ if one disregards the 
LDME. 

Of course, one could question our assumption 
that <C>3 Sl ( 3 5 j 1] )> = ^ (27+1) |R(0)| 2 and rather fit 
{(hs^S™)). In both the 7/i/r and TGS) cases, the corre¬ 
sponding LDMEs would then approximately be 100 times 
smaller. In particular, the singlet transition would be 100 
times less probable that what one expects from the leptonic 
decay. This would be an unlikely and dramatic violation of 
factorisation which should have implications elsewhere. In 
particular, a pair produced at short distances with the same 
quantum number as the physical state, among these the 
colour, would have a much larger probability to be broken 
up before eventually hadronising than expected. 

Although it is not as obvious as in the NRQCD formu¬ 
lae of Eq. 10, where the hypotheses of the CEM are trans¬ 
lated into direct relations between CO and CS transition 
probabilities, the same should happen in Eq. 8 where all the 
colour configurations are summed over and then considered 
on the same footage. In a process where the CS configura¬ 
tions dominate, such as qq —» y* —> QQ, CSM and CEM 
predictions necessarily differ. Contrary to NRQCD which 
encompasses the CSM, the CEM does not encompass the 
CSM. If one agrees with the data, the other cannot. The mat¬ 
ter is then how precise the predictions and the data are to rule 
out one approach or the other. 

Overall, one has to acknowledge that the conventional 
CEM central curves -as simplistic as the underlying idea of 
the model can be- give an account (Fig. 4 (a) & (b)) of the 
world data points as satisfactory as the LO CSM. The latter 
seems to underestimate the data at low energies while the 
former only has trouble to account for the TeV J/ijj points; 
the slope being more problematic than the normalisation 
which can be adjusted. All this is qualitative since the theo¬ 
retical uncertainties on the CEM are as large as that on the 
open heavy-flavour production which are admittedly large 
(see [27] for an up-to-date discussion on the cc case). 


7 Energy dependence of the colour-singlet channel at 
Born and one-loop accuracy 

As we just stated, the LO CSM curves are providing a sur¬ 
prisingly good description of the 7/i/f and T data at high 
energies without adjusting -and even less twisting- any pa¬ 
rameters. Although the central LO CSM curves agree with 
the data, the conventional theoretical uncertainties -from 
the arbitrary scales and the heavy-quark mass- are large 
(see the lower panels of Fig. 2)). It is therefore very natu¬ 
ral to look whether these uncertainties are reduced at one- 
loop accuracy. Such an observation was already made for 
the T case in [28] but this study was limited to a single V s , 
i.e. 200 GeV. 


7.1 Spin-triplet quarkonia: J/if/ and T 

Contrary the CO channels, the one-loop corrections to the 
CS channels only arise at (see e.g. Fig. lj & lk). Nev¬ 
ertheless, these are know since 2007 [76] and can also be 
computed with FDC as done in [77], In particular, there is 
no specific difficulty to integrate the a 3 and contributions 
in Pj since they are finite at Pr = 0. 

However, as already noted in [29], such NLO results 
tends to shows negative values at low Pj which can have a 
non-negligible impact on the total (i.e. Pj -integrated) cross 
section. To our knowledge, the energy dependence of the 
CSM at NLO has never been studied in detail. This is done 
below. 

Figs. 5 show the energy dependence of the NLO CSM (7 
curves). Note that, if a curve is not shown until 14 TeV, this 
indicates that the total yield got negative. The 3 red curves 
correspond to the default scale choices (^r =///■ = 2m q) and 
are indicative of the heavy-quark mass uncertainty, on the 
order of a factor of 4 for the 7/i/r and 2.5 for the T. In the 
former case, all 3 curves end up to be negative somewhere 
between 500 GeV and 2 TeV. Note also that the upper curve 
at low energies, i.e. for m c = 1.4 GeV, is the first to get nega¬ 
tive and crosses the other ones as if the negative contribution 
were more important for lighter systems 13 . In the T case, 
these 3 curves do not become negative at high energies -we 
have checked it up to ffs - 100 TeV. Nonetheless, they start 
to significantly differ from the LO curves (3 blue curves) 
above 1 TeV, contrary to the good LO vs NLO convergence 
found at RHIC energies in [28], One might thus be tempted 
to identify this weird energy behaviour to a low-x effect. 

Going further in the 7/i/r case, one can vary the factori¬ 
sation and renormalisation scales about the default choice. 
Doing so, one obtains two classes of curves. For /jr > fif 
(pink and orange), the yield remains positive, but it is not 
less unphysical for it to be practically constant as the energy 

13 unless the origin of this effect is due to gf > 2iiiq, see section 7.2. 
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Fig. 5: The cross section for direct (a) J/ip and (b) TILS’) as 
a function of cms-energy in the CSM at LO and NLO for 
various choices of the mass and scales compared with the 
existing experimental measurements (see text). 


increases between 1 and 10 TeV ! The only way to recover 
a semblance of increase is to take a large value of pr -and 
seemingly also a small value for pp. Obviously, whatever 
the reason for this behaviour is, for large enough pr, the 
QCD corrections which are proportional to a s (pR) neces¬ 
sarily get smaller and any difference between LO and NLO 
results should decrease. In the T case, as for the J/i// case, 
both curves with /jr > pp (pink and orange) correspond to 
the highest yields at high energies and the lowest at low en¬ 
ergies. When one chooses pp > /jr (purple and green), the 
high-energy yields become negative both for J/ip and T. In 
many respect, these observations are very similar to those 
made on the 1 5®. Such a pathological behaviour may thus 
not be related to the nature to the final state (see also next 
section). 

Large NNLO corrections are expected to show up at 
large energies (low x) as discussed in [78]. It is not clear 
if they could provide a solution to this issue. Another way to 


solve this might be to resum initial state radiation as done in 
the CEM [79] and for some CO channels [16]. 

At the light of such results, the most that one can reason¬ 
ably say is that the NLO CSM results may be reliable for T 
up to 200 GeV and for J/ip up to 60 GeV, that is up to \fs 
about 20 times the quarkonium mass. Above these value, the 
best that we have is the Born order results. 


7.2 Spin-singlet pseudo-scalar quarkonia: q c and t]b 

Contrary to the spin-triplet case, one can obtain analytical 
formulae [22, 80] for the spin-singlet pseudo-scalar produc¬ 
tion cross section such as that of rj c and This can in 
principle be of some help to understand the weird energy 
behaviour of the CS '.S’ \ yield and of some CO channels. 
Indeed, the LO production occurs as for some CO channels 
without final-state-gluon radiation. In fact the final state is 
simply colourless. 

As can be seen on Figs. 6, the issue is similar in many 
respects but for the fact that one does not obtain negative 
yields for the rjb- For the rj c , the curves for /jr > pp remains 
positive at high energies -as for the 7/ ip. One also sees that 
the crossing of the central LO and NLO curves occurs at 
larger \[s than for the ’’.S' i states. However, such small quan¬ 
titative differences may be due to the fact that we computed 
the y-integrated cross sections using the analytical expres¬ 
sions of [22] instead of the y-differential cross section at 

>’ = o. 

We have investigated this in more details by looking at 
the different NLO contributions (the real emissions from gg 
and qg fusion as well as the virtual (loop) contributions) in 
order to see which channels induce the negative contribu¬ 
tions and for which scale/mass values. However, it must be 
stressed that the decomposition between these different con¬ 
tributions depend on the regularisation method used. For in¬ 
stance, the decomposition is drastically different when us¬ 
ing FDC -with sometimes a very large cancellation between 
the positive real-emission gg contribution and the negative 
sum of the Born and loop gg contributions - and the formu¬ 
lae of [22, 80] -where all the gg contributions are gathered. 
Yet, we checked that we obtain exactly the same results with 
both methods; the regularisation method or numerical insta¬ 
bilities cannot be the source of the issues observed above. 

As regards the qg contribution, only pp/niQ matters to 
tell whether it will change sign. For pp close to ihq and be¬ 
low, it will be positive (negative) at small (large) sfs. Other¬ 
wise, it remains negative for any sfs. The value of pr only 
influences the normalisation. 

As what concerns the gg contributions, which are ex¬ 
pected to be dominant at high energies, both pp/niQ 
and pr/pf matter. For pp - ihq, the gg contribution 
monotonously increase as function sfs irrespective of 
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Fig. 6 : The cross section for direct (a) q r and (b) ?//, as a func¬ 
tion of cms-energy in the CSM at LO and NLO for various 
choice of the mass and scales. 


A'r/A'f- For /_//. > mg , the gq one is rather small and, despite 
being negative, does not induce a turn over in the increase 
of the cross section. For - 2ihq, the gg contribution gets 
negative at large Vs for /jr < qip. For /Ur > p/, it remains al¬ 
ways positive. Yet the sum gg+gq can still become negative 
since, in some cases, gq increases faster with \fs. 

All this seems in lines, for the gg [gq] part, with the for¬ 
mulae (C.25) and (C.26) [(C.32) and (C.35)] in the appendix 
C of [22]. Both contributions indeed exhibit logarithms of 
/jp/mQ multiplied by a factor function of ihq/s. 

As we can see, the results are already difficult to inter¬ 
pret for the spin singlet case. For the spin triplet, we do not 
have similar analytical results. We can only guess that the 
stmcture is similar. 

What seems surprising is that when one inspects simi¬ 
lar expressions for the i]q production at NLO in the TMD 
factorisation [32], such negative terms do not appear as ob¬ 


vious. We are thus entitled to wonder whether such a for¬ 
malism, which automatically resums logarithm of transverse 
momenta, may provide the solution to this issue. Another 
possible solution may be the consideration of NNLO cor¬ 
rections, which may show opposite signs to that at NLO. 
This is obviously beyond the scope of our analysis. 


8 Conclusion 

We have performed the first full analysis of the energy 
dependence of the quarkonium-production cross section at 
one-loop accuracy both in NRQCD and in the CSM. Taken 
at face value, our results would indicate a severe break down 
of NRQCD universality - in line with the previous analysis 
of Maltoni et al. [15] - unless one keeps the LDMEs close 
to the fit of Butenschoen and Kniehl, which however dis¬ 
agrees with the J/if/ polarisation measurements and the q c 
cross sections. 

The situation is however slightly more intricate since we 
have uncovered a weird -sometimes unphysical- behaviour 
at large energies where one approaches the small-jc regime 
where non-linear effects in the parton densities may be rel¬ 
evant. This certainly casts doubts on the numerical values 
we obtained at LHC energies since collinear factorisation 
on which we based our analysis could break down. 

Yet, up to a few hundred GeV, the energy dependence of 
the different octet channels at NLO seems well behaved and 
there is no reason to doubt this. In this region, the NLO yield 
prediction by NRQCD after fitting the mid- and h igh-/ J y 
quarkonium data -i.e. the yield and its polarisation- would 
overshoot the data by a factor ranging from 4 to more than 
10. The same holds true at LO (see Appendix A). To repro¬ 
duce the data, the CO LDMEs should be much smaller than 
what they are found in order to be to reproduce the Tevatron 
and LHC P 7 -differential cross section in the case of the J/ip 
and i[s(2S ). 

On the other hand, the LO CSM provides a decent en¬ 
ergy dependence in agreement with the existing data, ex¬ 
cept for V?<40 GeV and is therefore up to ten times be¬ 
low the -data-overshooting- CO contributions. At one loop, 
the results are however ill-behaved for the charmonia and 
the cross sections can even become negative at large -\fs for 
some -reasonable- scale choices. The situation is a bit better 
for T(l.S'). The same occurs for the spin-singlet quarkonia. 
In this case, one-loop results exist in the TMD factorisation 
approach and do not seem to be prone to such an issue. In 
the case of double J /<// production, the energy dependence 
at one loop seems well-behaved [81], Finally, we investi¬ 
gated the energy dependence of the yield in the CEM where 
the final states are treated rather differently and we did not 
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find any specific problems 14 . We are therefore tempted to 
attribute this problem to initial-state effects. 

In [82], Ma and Venugopalan obtained a good descrip¬ 
tion of the low-/ J y //(A data over a wide range of energy by, 
on the one hand, using the LDMEs from [17] -our second 
set- and, on the other, a CGC-based computation of the low- 
Pj dependence. In reproducing the data, they found that the 
CS contribution is only 10% of the total yield. This 10% is 
reminiscent of the factor 10 between the CS and CO in our 
“collinear” study. From our viewpoint, it looks as if the spe¬ 
cific ingredient of this CGC-based computation would cor¬ 
respond to an effective reduction of the two-gluon flux 15 by 
a factor of 10. It is therefore very interesting to find out new 
processes which would be sensitive to this physics. 

The negative yields obtained in the collinear case - ob¬ 
served for 1 5 1 , 'Sq'V 3 S [ / ] , 3 / j ^ channels- could also be 

cured by adding the large contributions of the one-loop am¬ 
plitude squared -thus positive. This may look like an ad- 
hoc solution which certainly questions the convergence of 
the perturbative series in a s . However, large NNLO correc¬ 
tions have already been discussed 10 years ago in [78]. An¬ 
other path towards a solution may also be higher-twist con¬ 
tributions [83] where two gluons come from a single pro¬ 
ton as recently rediscussed in [84]. In any case, whatever 
the explanation for this situation may be, past claims that 
colour-octet transitions are dominantly responsible for low- 
Pj quarkonium production were premature in light of the 
results presented here. 
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Appendix A: LO NRQCD predictions 

Sharma and Vitev recently performed [85] a LO fit of the 
CO LDMEs using RHIC, Tevatron and LHC J/i// data. Set¬ 
ting (Oj/f( l S l * ] )) = (Oj/^P^/ml and accounting for the 
possible feed-downs, they obtained: 

- q 81 )) = 0.018 GeV 3 , 

- <<9^( 3 5 [ 1 8] )) = 0.0013 GeV 3 . 

14 Apart from the fact that the recasting of the CEM into NRQCD does 
not seem to work, phenomenologically speaking. 

15 The comparison is however a bit more complex since this CGC- 
based approach accounts for contributions which are normally sup¬ 
pressed in the collinear limit. 
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Fig. 7: The cross section for the production of a J/iJj at LO 
from only colour octet states as a function of the cms-energy 
for various choice of the mass and scales. 


Based on this fit, we derive the energy dependence for 
the direct J/i// which is shown on Fig. 7. Without any sur¬ 
prise, the results badly overshoot the world data. 
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